#0) Targets (uses mix.var and U0 from environment)
target <- function(num.v){
olydum <- n.firms<n.models
u0dum <- U0==1
ogpct <- fifelse(OGpct == 50,50,91)
DT <- data.table(vnum =1:5,
                       eta.target = rep(4,5),
                       price.elas.target = rep(-4,5),
                       ISE.target = c(0,0,NA,0.85,0.85),
                       CR5.target = rep(86,5), # BLP data for 1990 (CR5 of firms)
                       dom.share.target = rep(68,5), # BLP data for 1990 (dom.share)
                       out.share.target = rep(ogpct,5)*u0dum + rep(0,5)*(1-u0dum),
                       crosselasmaha.target = c(0,-0.53,-0.97,-0.97,-0.97),
#                      meanmaha.target = c(2.53,2.53,2.81,2.81,2.81),
#                      IQR.target = c(1.79,1.79,1.97,1.97,1.79),
                       exvar.target = c(0,0.01,0.017,0.017,0.01),
                       mncv.target = c(0,0.523,0.885,0.885,0.523)
)
if(olydum==0) DT[, CR5.target := rep(18,5)] # BLP data for 1990 (CR5 of models)
return(DT[vnum <= num.v])
}
#
#1) DGP for 1 rep
#eqbm for one market
eqbm.calib <- function(settings,v=1,market=1) { 
  #exogenous part (needs to also have trade costs matrix because of dave)
  tau.grid.allv <-  as.data.table(make.tradecosts(n.countries,n.firms,n.models,v,settings=settings))
  tau.grid.allv[,tariff.change:=-tar] # new tar.change is a dummy, tariff.change is a variable with desired changes.
  data.exog.mm <<- DGP.exog(settings=settings,v=v,grid=subset(tau.grid.allv,vtau==v,select=-vtau)) #gen exog data
  #equilibrium + CES approx 
  if(substr(mix.var,1,4)=="MCES"){
  eqbm.pre.mm  <- as.data.table(DGP.endog(parallel=F,v=v,tar.change=0,data.exog = data.exog.mm)) #eqbm
  } else eqbm.pre.mm  <- as.data.table(DGP.endog.mlog(parallel=F,v=v,tar.change=0,data.exog = data.exog.mm))
  eqbm.pre.mm.with.ces <- as.data.table(estapprox.mm(eqbm.pre.mm)) #add ces approximation results
  #oligopoly version 
  if(substr(mix.var,6,8)=="OLY"){
    eqbm.pre.mm.with.ces[, s.f := sum(s), by=c("f","n")]
    eta.ces.oly.init <<- 3
    FPIresults <- fpiter(eta.ces.oly.init,fixptfn=ces.oly.mm.update,pre=eqbm.pre.mm.with.ces,nu=1) #estimate eta.ces.oly
    eta.ces.oly <- FPIresults$par 
    eqbm.pre.mm.with.ces[, eta := eta.ces.oly]
  }
  pre <- eqbm.pre.mm.with.ces[n==market,] #keep only one destination
  return(pre)
  }
#collecting moments 
  monte.func.mm.calib <- function(repi,settings,v=1,market=1,cf.meth="EHA") { # monte carlo stage
  #
  #compute pre equilibrium
  pre <- eqbm.calib(settings,v=v,market=market)
  #
  #richsubs moments 
  EQ <- as.data.frame(pre)
  if(v==3 | v==4)  {
    MD <- MahaDist(data=EQ,features = c(paste0("x",1:n.x),"p")) # including p should be linked to the version
  } else MD <- MahaDist(data=EQ,features = c(paste0("x",1:n.x)))
  if(substr(mix.var,1,4)=="MCES"){
  XE <- rich.subs(dat=data.exog.mm,eqbm.data = EQ,market=market) 
  } else   XE <- rich.subs.mlog(dat=data.exog.mm,eqbm.data = EQ,market=market) 
  DT <- merge(XE,MD,by=c("m","j"))
  res_loglin <- feols(log(cross_elas)~ maha_dist | m+j,DT[m !=j])  
  crosselasmaha <- res_loglin$coefficients[[1]] 
  exvar <- DT[m==j, mean(wvp_m/s_m)] 
  DTcross <- DT[m!=j]
  DTcross <- DTcross[,.(sd.cross_elas_m=sd(cross_elas),mn.cross_elas_m=mean(cross_elas)), by=.(m)]
  mncv <- DTcross[,mean(sd.cross_elas_m/mn.cross_elas_m)] 
  #
  #other moments
  out.shr <- with(pre,1-sum(s)) # initial outside good share 
  pre$s <- pre$s/sum(pre$s)  # those should be only on INSIDERS
  cr5 <- CR5.agg(pre)
  dom.shr <- with(pre,sum(s[i==n])) # initial domestic share 
  price.elas <- with(pre,mean(pelas.m))  # mean (across models) of own price elasticity
  #convergence diagnostic
  converge.prob <- max(pre$n.loops)==MAXIT 
  #return
  return(data.frame(repi=repi,v=v, eta=pre$eta[1], price.elas, ISE=pre$ISE[1], cr5, dom.shr, out.shr, crosselasmaha, exvar, mncv, converge.prob))
}
#
# 2) Doing multiple monte carlo reps for a single market
monte.rep.mm.calib <- function(parallel=TRUE,settings,v=1,market=1,cf.meth="EHA"){ #get equilibrium market shares and prices for each market and append them
  if(parallel==T)  foreach(i=1:n.reps, .combine = 'rbind')  %dorng%  monte.func.mm.calib(i,settings,v,market,cf.meth) 
  else foreach(i=1:n.reps, .combine = 'rbind')  %do%  monte.func.mm.calib(i,settings,v,market,cf.meth) 
}
#
# 3) Moments (rounded to the desired level of tolerance)
monte.moments.mm <- function(simdata,round.moments=TRUE){ # simdata is a data.table
  simdata <-simdata[, converge.prob.mean := mean(converge.prob),by=.(v)]
 if(round.moments==TRUE) {
   monte.summary <-simdata[converge.prob==FALSE,.(eta=round(mean(eta),digits=1), #all
                              price.elas=round(mean(price.elas),digits=1),
                              ISE=round(mean(ISE),digits=2), # only for settings 4 and 5
                             cr5=round(100*mean(cr5),digits=0), # panel b and c only but all settings
                             dom.shr=round(100*mean(dom.shr),digits=0),
                             out.shr=round(100*mean(out.shr),digits=0),
                             crosselasmaha=round(mean(crosselasmaha),digits=2),
                             exvar=round(mean(exvar),digits=3),
                             mncv=round(mean(mncv),digits=3),
                             converge.prob.mean=mean(converge.prob.mean)) #all
                          ,by=.(v)
                          ]
 } else {
   monte.summary <-simdata[converge.prob==FALSE,.(eta=mean(eta), #all
                              price.elas=mean(price.elas),
                              ISE=mean(ISE), # only for settings 4 and 5
                              cr5=100*mean(cr5), # panel b and c only but all settings
                              dom.shr=100*mean(dom.shr),
                              out.shr=100*mean(out.shr),
                              crosselasmaha=mean(crosselasmaha),
                              exvar=mean(exvar),
                              mncv=mean(mncv),
                              converge.prob.mean=mean(converge.prob.mean)) #all
                           ,by=.(v)
                           ]
 }
  #
  if(sum(monte.summary$converge.prob.mean)>0) cat("Alert: at least one replication hit max iterations")
  #
  return(monte.summary) # data.table with 5 rows (1 for each v), 5 columns (1 for each targeted moment)
}

# 3') Outputs the moment
monte.moments.output <- function(monte.moments.data){
monte.mom <- copy(monte.moments.data)
sink(file=paste("../Tables/CF_MM/",mix.var,"/CF_calib_",pan.var,"OG_",OGpct,".txt",sep=""))
cat(paste(mix.var,pan.var,": Number of reps ",n.reps,"\n","\n"))
cat(paste("Moments: ","\n"))
print( monte.mom )
monte.mom[,c("converge.prob.mean") := NULL]
cat(paste("\n","\n"))
cat(paste("Moments - Targets: ","\n"))
print( monte.mom - target(num.v=nv) )
sink()
}



